Conjugate Sampler - Split/Merge Demo
Jonathan Klus
2023-12-7
DEE Sampler Setup
- \(y_1, \dots, y_n \sim \sum_{k=1}^K\pi_kN(\mu_k, \sigma^2 I_J)\)
- \(\mu_1, \dots, \mu_K \sim N(\mu_0, r\sigma^2 I_J)\)
- \(\sigma^2 \sim IG(a,b)\)
- \(r \sim IG(g, h)\)
- \(b \sim Ga(d, f)\)
DEV Sampler Setup
- \(y_1, \dots, y_n \sim \sum_{k=1}^K\pi_kN(\mu_k, \sigma^2_k I_J)\)
- \(\mu_1, \dots, \mu_K \sim N(\mu_0, r\sigma^2_k I_J)\)
- \(\sigma^2_1,\dots, \sigma^2_K \sim IG(a,b)\)
- \(r \sim IG(g, h)\)
- \(b \sim Ga(d, f)\)
Sampler Examples
2D Example - well separated
Mixture of 3 multivariate normal densities all with diagonal variance structure and \(\sigma^2 = 10\). The true means are \((-20,20)\), \((20,-20)\), and \((0,0)\), and \(\sum_{j}n_j = 30\).
DEV sampler
without split-merge step
# n_iter = 5000
# use_cores = min(parallel::detectCores(), nreps)
# sim_results_DEV_no_sm = parallel::mclapply(X = 1:10, mc.cores = use_cores,
# FUN = function(x){MVN_CRP_sampler_DEV(
# S = n_iter, seed = seeds[x], y = yreps[[x]]$y,
# alpha = 1, r = 10, g = 1, h = 50,
# sigma_hyperprior = FALSE, fix_r = FALSE,
# mu0 = matrix(round((colMeans(matrix(unlist(yreps[[x]]$y), ncol = 2))),0), ncol = 1),
# a = 1, b = 50,
# k_init = 1, diag_weights = FALSE,
# verbose = FALSE, split_merge = FALSE
# )}
# )
#
# saveRDS(object = sim_results_DEV_no_sm, file = "../MCMC_Runs/conjDEVsamp_minisimstudy_noSM_6_Dec_2023.rds")
sim_results_DEV_no_sm = readRDS(file = "../MCMC_Runs/conjDEVsamp_minisimstudy_noSM_6_Dec_2023.rds")
for(rep in 1:nreps){
print(table(sim_results_DEV_no_sm[[rep]]$k))
p1 = make_k_traceplot(k = sim_results_DEV_no_sm[[rep]]$k,
group_assign = sim_results_DEV_no_sm[[rep]]$group_assign)
print(p1)
}##
## 1 2 3 4 5 6 7
## 2 15 3965 856 131 24 7
##
## 1 2 3 4 5
## 4 7 4535 428 26
##
## 1 3 4 5 6 7 8
## 1 1694 2980 295 24 5 1
##
## 1 3 4 5 6 7
## 2 4447 478 61 10 2
##
## 1 2 3 4 5 6 7
## 5 1 459 3600 805 119 11
##
## 1 2 3 4 5 6 7
## 6 2 1632 3008 335 11 6
##
## 1 2 3 4 5 6 7
## 1 4 4094 818 70 11 2
##
## 1 2 3 4 5 6 7
## 4 4 2260 2310 373 44 5
##
## 1 2 3 4 5 6
## 4 5 3833 1079 76 3
##
## 1 2 3 4 5 6
## 5 44 3902 924 111 14
with split-merge step
# sim_results_DEV_with_sm = parallel::mclapply(X = 1:10, mc.cores = use_cores,
# FUN = function(x){MVN_CRP_sampler_DEV(
# S = n_iter, seed = seeds[x], y = yreps[[x]]$y,
# alpha = 1, r = 10, g = 1, h = 50,
# sigma_hyperprior = FALSE, fix_r = FALSE,
# mu0 = matrix(round((colMeans(matrix(unlist(yreps[[x]]$y), ncol = 2))),0), ncol = 1),
# a = 1, b = 50,
# k_init = 1, diag_weights = FALSE,
# verbose = FALSE, split_merge = TRUE, sm_iter = 10
# )}
# )
#
# saveRDS(object = sim_results_DEV_with_sm, file = "../MCMC_Runs/conjDEVsamp_minisimstudy_withSM_6_Dec_2023.rds")
sim_results_DEV_with_sm = readRDS(file = "../MCMC_Runs/conjDEVsamp_minisimstudy_withSM_6_Dec_2023.rds")
for(rep in 1:nreps){
print(table(sim_results_DEV_with_sm[[rep]]$k))
p1 = make_k_traceplot(k = sim_results_DEV_with_sm[[rep]]$k,
group_assign = sim_results_DEV_with_sm[[rep]]$group_assign)
print(p1)
}##
## 1 2 3 4 5 6 7
## 2 5 3399 1243 309 41 1
##
## 1 2 3 4 5 6
## 5 11 4512 440 29 3
##
## 1 2 3 4 5 6 7
## 1 3 4302 603 78 10 3
##
## 1 2 3 4 5 6 7
## 2 55 3601 1117 196 26 3
##
## 1 2 3 4 5 6 7
## 5 1 3743 1108 127 15 1
##
## 1 2 3 4 5 6 7 8
## 6 2 3504 1207 248 27 5 1
##
## 1 2 3 4 5 6
## 1 4 3627 1205 152 11
##
## 1 2 3 4 5 6
## 4 4 3288 1521 158 25
##
## 1 3 4 5 7
## 4 4691 299 5 1
##
## 1 2 3 4 5 6 7
## 5 25 4074 765 120 10 1
DEE sampler
without split-merge step
# sim_results_DEE_no_sm = parallel::mclapply(X = 1:10, mc.cores = use_cores,
# FUN = function(x){MVN_CRP_sampler_DEE(
# S = n_iter, seed = seeds[x], y = yreps[[x]]$y,
# alpha = 1, r = 10, g = 1, h = 50,
# sigma_hyperprior = FALSE, fix_r = FALSE,
# mu0 = matrix(round((colMeans(matrix(unlist(yreps[[x]]$y), ncol = 2))),0), ncol = 1),
# a = 1, b = 50,
# k_init = 1, diag_weights = FALSE,
# verbose = FALSE, split_merge = FALSE
# )}
# )
#
# saveRDS(object = sim_results_DEE_no_sm, file = "../MCMC_Runs/conjDEEsamp_minisimstudy_noSM_6_Dec_2023.rds")
sim_results_DEE_no_sm = readRDS(file = "../MCMC_Runs/conjDEEsamp_minisimstudy_noSM_6_Dec_2023.rds")
for(rep in 1:nreps){
print(table(sim_results_DEE_no_sm[[rep]]$k))
p1 = make_k_traceplot(k = sim_results_DEE_no_sm[[rep]]$k,
group_assign = sim_results_DEE_no_sm[[rep]]$group_assign)
print(p1)
}##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 1 1 7 23 94 230 450 669 710 712 610 496 351 280 172 111 49 17 14 3
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 1 2 14 71 216 449 620 740 715 682 518 333 262 176 115 50 20 12 2 1
## 21
## 1
##
## 1 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
## 1 11 67 261 538 724 739 666 593 392 301 234 170 138 80 46 26 8 4 1
##
## 1 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 2 1 13 91 196 359 574 686 756 650 516 457 311 191 101 66 21 7 2
##
## 1 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 1 6 43 144 327 523 627 734 640 611 494 331 261 131 67 41 14 4 1
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 3 3 1 20 92 277 515 723 805 680 626 419 325 206 147 90 49 12 6 1
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 1 1 5 42 146 387 620 845 800 722 488 338 245 158 102 56 24 14 3 1
## 21
## 2
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
## 4 3 2 25 126 376 621 804 787 698 522 387 280 159 95 49 36 22 4
##
## 1 3 4 5 6 7 8 9 10 11 12 13 14 15 16
## 4 61 217 445 677 780 781 728 538 333 207 118 77 22 12
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 5 3 9 32 109 259 468 680 765 715 574 476 347 231 135 85 53 39 10 5
with split-merge step
# sim_results_DEE_with_sm = parallel::mclapply(X = 1:10, mc.cores = use_cores,
# FUN = function(x){MVN_CRP_sampler_DEE(
# S = n_iter, seed = seeds[x], y = yreps[[x]]$y,
# alpha = 1, r = 10, g = 1, h = 50,
# sigma_hyperprior = FALSE, fix_r = FALSE,
# mu0 = matrix(round((colMeans(matrix(unlist(yreps[[x]]$y), ncol = 2))),0), ncol = 1),
# a = 1, b = 50,
# k_init = 1, diag_weights = FALSE,
# verbose = FALSE, split_merge = TRUE, sm_iter = 10
# )}
# )
#
# saveRDS(object = sim_results_DEE_with_sm, file = "../MCMC_Runs/conjDEEsamp_minisimstudy_withSM_6_Dec_2023.rds")
sim_results_DEE_with_sm = readRDS(file = "../MCMC_Runs/conjDEEsamp_minisimstudy_withSM_6_Dec_2023.rds")
for(rep in 1:nreps){
print(table(sim_results_DEE_with_sm[[rep]]$k))
p1 = make_k_traceplot(k = sim_results_DEE_with_sm[[rep]]$k,
group_assign = sim_results_DEE_with_sm[[rep]]$group_assign)
print(p1)
}##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 1 1 14 30 87 273 468 725 750 663 580 492 365 221 161 96 47 17 6 3
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 1 2 78 111 259 432 632 690 748 652 441 351 281 135 81 66 23 8 7 2
##
## 1 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
## 1 24 69 250 536 781 788 672 539 406 272 231 161 118 67 46 22 12 3 2
##
## 1 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
## 2 4 41 164 371 600 755 706 660 494 410 275 226 141 84 38 19 6 1 2
## 22
## 1
##
## 1 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 1 4 39 151 336 541 682 689 675 553 458 358 222 144 83 32 18 12 2
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 3 3 4 26 114 300 541 716 820 688 560 425 309 209 136 76 45 17 7 1
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 1 1 4 29 180 403 718 855 755 627 498 320 236 167 84 67 35 13 3 4
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 4 3 7 48 189 413 647 805 762 673 500 355 244 149 99 55 25 13 8 1
##
## 1 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
## 4 83 192 494 679 762 828 654 522 329 197 131 78 33 9 4 1
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 5 3 4 43 137 345 556 700 741 658 514 430 303 201 148 96 62 30 19 5
UVV sampler
without split-merge step
sim_results_UVV_no_sm = parallel::mclapply(X = 1:10, mc.cores = use_cores,
FUN = function(x){MVN_CRP_sampler_UVV(
S = n_iter, seed = seeds[x], y = yreps[[x]]$y,
alpha = 1, r = 10, g = 1, h = 50, nu = 2,
fix_r = FALSE,
mu0 = matrix(round((colMeans(matrix(unlist(yreps[[x]]$y), ncol = 2))),0), ncol = 1),
lambda0 = diag(x = 15, nrow = 2),
k_init = 1, diag_weights = FALSE,
verbose = FALSE, split_merge = FALSE
)}
)
saveRDS(object = sim_results_UVV_no_sm, file = "../MCMC_Runs/conjUVVsamp_minisimstudy_noSM_6_Dec_2023.rds")
# for(rep in 1:nreps){
# print(table(sim_results_DEE_with_sm[[rep]]$k))
# p1 = make_k_traceplot(k = sim_results_DEE_with_sm[[rep]]$k,
# group_assign = sim_results_DEE_with_sm[[rep]]$group_assign)
# print(p1)
# }with split-merge step
sim_results_UVV_with_sm = parallel::mclapply(X = 1:10, mc.cores = use_cores,
FUN = function(x){MVN_CRP_sampler_UVV(
S = n_iter, seed = seeds[x], y = yreps[[x]]$y,
alpha = 1, r = 10, g = 1, h = 50, nu = 2,
fix_r = FALSE,
mu0 = matrix(round((colMeans(matrix(unlist(yreps[[x]]$y), ncol = 2))),0), ncol = 1),
lambda0 = diag(x = 15, nrow = 2),
k_init = 1, diag_weights = FALSE,
verbose = FALSE, split_merge = TRUE, sm_iter = 10
)}
)
saveRDS(object = sim_results_UVV_with_sm, file = "../MCMC_Runs/conjUVVsamp_minisimstudy_withSM_6_Dec_2023.rds")2D Example - close together
Mixture of 3 multivariate normal densities all with diagonal variance structure and \(\sigma^2 = 10\). The true means are \((10,10)\), \((0,-5)\), , \((12,-2)\), and \((0,10)\), and \(\sum_{j}n_j = 50\).
DEV sampler
without split-merge step
# n_iter = 5000
# use_cores = min(parallel::detectCores(), nreps)
# sim_results_DEV_no_sm = parallel::mclapply(X = 1:10, mc.cores = use_cores,
# FUN = function(x){MVN_CRP_sampler_DEV(
# S = n_iter, seed = seeds[x], y = yreps[[x]]$y,
# alpha = 1, r = 10, g = 1, h = 50,
# sigma_hyperprior = FALSE, fix_r = FALSE,
# mu0 = matrix(round((colMeans(matrix(unlist(yreps[[x]]$y), ncol = 2))),0), ncol = 1),
# a = 1, b = 50,
# k_init = 1, diag_weights = FALSE,
# verbose = FALSE, split_merge = FALSE
# )}
# )
#
# saveRDS(object = sim_results_DEV_no_sm, file = "../MCMC_Runs/conjDEVsamp_minisimstudy_close_noSM_6_Dec_2023.rds")
sim_results_DEV_no_sm = readRDS(file = "../MCMC_Runs/conjDEVsamp_minisimstudy_close_noSM_6_Dec_2023.rds")
for(rep in 1:nreps){
print(table(sim_results_DEV_no_sm[[rep]]$k))
p1 = make_k_traceplot(k = sim_results_DEV_no_sm[[rep]]$k,
group_assign = sim_results_DEV_no_sm[[rep]]$group_assign)
print(p1)
}##
## 1 2 3 4 5 6 7 8 9
## 11 317 2558 1416 522 143 26 6 1
##
## 1 2 3 4 5 6 7 8
## 62 469 2370 1431 541 109 15 3
##
## 1 2 3 4 5 6 7 8 9
## 250 104 3296 1033 228 67 19 2 1
##
## 1 2 3 4 5 6 7 8
## 123 123 2814 1295 504 120 19 2
##
## 1 2 3 4 5 6 7 8 9
## 230 114 1643 1782 851 302 66 11 1
##
## 1 2 3 4 5 6 7 8
## 111 150 3301 1154 237 40 5 2
##
## 1 2 3 4 5 6 7 8
## 99 60 3836 795 152 50 7 1
##
## 1 2 3 4 5 6 7 8 9
## 359 604 1080 1566 897 375 100 18 1
##
## 1 2 3 4 5 6
## 51 20 4043 783 100 3
##
## 1 2 3 4 5 6 7 8
## 202 522 2331 1254 496 149 40 6
with split-merge step
# sim_results_DEV_with_sm = parallel::mclapply(X = 1:10, mc.cores = use_cores,
# FUN = function(x){MVN_CRP_sampler_DEV(
# S = n_iter, seed = seeds[x], y = yreps[[x]]$y,
# alpha = 1, r = 10, g = 1, h = 50,
# sigma_hyperprior = FALSE, fix_r = FALSE,
# mu0 = matrix(round((colMeans(matrix(unlist(yreps[[x]]$y), ncol = 2))),0), ncol = 1),
# a = 1, b = 50,
# k_init = 1, diag_weights = FALSE,
# verbose = FALSE, split_merge = TRUE, sm_iter = 10
# )}
# )
#
# saveRDS(object = sim_results_DEV_with_sm,
# file = "../MCMC_Runs/conjDEVsamp_minisimstudy_close_withSM_6_Dec_2023.rds")
sim_results_DEV_with_sm = readRDS(file = "../MCMC_Runs/conjDEVsamp_minisimstudy_close_withSM_6_Dec_2023.rds")
for(rep in 1:nreps){
print(table(sim_results_DEV_with_sm[[rep]]$k))
p1 = make_k_traceplot(k = sim_results_DEV_with_sm[[rep]]$k,
group_assign = sim_results_DEV_with_sm[[rep]]$group_assign)
print(p1)
}##
## 1 2 3 4 5 6 7 8
## 1370 998 1134 920 401 145 26 6
##
## 1 2 3 4 5 6 7
## 2388 1425 795 275 95 19 3
##
## 1 2 3 4 5 6 7
## 2442 1132 789 457 147 30 3
##
## 1 2 3 4 5 6 7 8
## 1860 1193 1110 569 189 62 15 2
##
## 1 2 3 4 5 6 7
## 2857 1049 667 279 123 19 6
##
## 1 2 3 4 5 6 7 8
## 1515 1005 1300 730 320 105 19 6
##
## 1 2 3 4 5 6 7
## 3004 1250 502 174 49 18 3
##
## 1 2 3 4 5 6 7 8
## 2751 1083 629 355 125 40 14 3
##
## 1 2 3 4 5 6 8
## 1946 1168 1356 429 94 6 1
##
## 1 2 3 4 5 6 7
## 2193 1296 876 405 170 51 9
DEE sampler
without split-merge step
# sim_results_DEE_no_sm = parallel::mclapply(X = 1:10, mc.cores = use_cores,
# FUN = function(x){MVN_CRP_sampler_DEE(
# S = n_iter, seed = seeds[x], y = yreps[[x]]$y,
# alpha = 1, r = 10, g = 1, h = 50,
# sigma_hyperprior = FALSE, fix_r = FALSE,
# mu0 = matrix(round((colMeans(matrix(unlist(yreps[[x]]$y), ncol = 2))),0), ncol = 1),
# a = 1, b = 50,
# k_init = 1, diag_weights = FALSE,
# verbose = FALSE, split_merge = FALSE
# )}
# )
#
# saveRDS(object = sim_results_DEE_no_sm, file = "../MCMC_Runs/conjDEEsamp_minisimstudy_close_noSM_6_Dec_2023.rds")
sim_results_DEE_no_sm = readRDS(file = "../MCMC_Runs/conjDEEsamp_minisimstudy_close_noSM_6_Dec_2023.rds")
for(rep in 1:nreps){
print(table(sim_results_DEE_no_sm[[rep]]$k))
p1 = make_k_traceplot(k = sim_results_DEE_no_sm[[rep]]$k,
group_assign = sim_results_DEE_no_sm[[rep]]$group_assign)
print(p1)
}##
## 1 2 3 4 6 8 9 10 11 12 13 14 15 16 17 18
## 15 4 1 1 1 2 21 115 343 729 1064 1016 787 486 258 103
## 19 20 21 22 23
## 36 13 3 1 1
##
## 1 2 3 4 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
## 29 2 1 2 1 3 18 72 225 540 756 952 944 706 431 209 68 32 7 2
##
## 1 2 3 5 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
## 11 2 3 1 3 15 44 172 421 733 895 918 718 477 317 145 84 24 11 4
## 23
## 2
##
## 1 2 3 4 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
## 2 2 1 1 1 1 8 49 196 508 831 980 917 650 419 227 114 59 24 9
## 22
## 1
##
## 1 2 3 5 7 8 9 10 11 12 13 14 15 16 17 18
## 8 10 1 1 1 2 14 89 358 706 1006 1036 773 518 271 130
## 19 20 21 22
## 53 16 5 2
##
## 1 2 3 4 6 7 8 9 10 11 12 13 14 15 16 17
## 2 1 1 1 1 1 7 43 213 465 832 1002 924 708 443 221
## 18 19 20 21 22 23
## 93 32 3 5 1 1
##
## 1 2 3 4 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
## 2 2 1 1 2 11 42 174 409 692 830 914 758 511 318 181 89 40 16 4
## 22 23
## 2 1
##
## 1 2 3 4 5 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
## 1 3 2 1 1 1 1 22 72 253 556 852 921 887 633 439 210 105 35 5
##
## 1 2 3 5 6 7 8 9 10 11 12 13 14 15 16 17
## 20 1 2 2 15 53 189 467 814 1005 885 692 444 253 112 34
## 18
## 12
##
## 1 2 3 4 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
## 8 1 1 1 1 6 28 133 363 660 909 992 807 541 333 145 49 19 1 1
## 24
## 1
with split-merge step
# sim_results_DEE_with_sm = parallel::mclapply(X = 1:10, mc.cores = use_cores,
# FUN = function(x){MVN_CRP_sampler_DEE(
# S = n_iter, seed = seeds[x], y = yreps[[x]]$y,
# alpha = 1, r = 10, g = 1, h = 50,
# sigma_hyperprior = FALSE, fix_r = FALSE,
# mu0 = matrix(round((colMeans(matrix(unlist(yreps[[x]]$y), ncol = 2))),0), ncol = 1),
# a = 1, b = 50,
# k_init = 1, diag_weights = FALSE,
# verbose = FALSE, split_merge = TRUE, sm_iter = 10
# )}
# )
#
# saveRDS(object = sim_results_DEE_with_sm, file = "../MCMC_Runs/conjDEEsamp_minisimstudy_close_withSM_6_Dec_2023.rds")
sim_results_DEE_with_sm = readRDS(file = "../MCMC_Runs/conjDEEsamp_minisimstudy_close_withSM_6_Dec_2023.rds")
for(rep in 1:nreps){
print(table(sim_results_DEE_with_sm[[rep]]$k))
p1 = make_k_traceplot(k = sim_results_DEE_with_sm[[rep]]$k,
group_assign = sim_results_DEE_with_sm[[rep]]$group_assign)
print(p1)
}Some iterations failed.
UVV sampler
without split-merge step
# sim_results_UVV_no_sm = parallel::mclapply(X = 1:10, mc.cores = use_cores,
# FUN = function(x){MVN_CRP_sampler_UVV(
# S = n_iter, seed = seeds[x], y = yreps[[x]]$y,
# alpha = 1, r = 10, g = 1, h = 50, nu = 2,
# fix_r = FALSE,
# mu0 = matrix(round((colMeans(matrix(unlist(yreps[[x]]$y), ncol = 2))),0), ncol = 1),
# lambda0 = diag(x = 15, nrow = 2),
# k_init = 1, diag_weights = FALSE,
# verbose = FALSE, split_merge = FALSE
# )}
# )
#
# saveRDS(object = sim_results_UVV_no_sm, file = "../MCMC_Runs/conjUVVsamp_minisimstudy_close_noSM_6_Dec_2023.rds")
# sim_results_UVV_no_sm = readRDS(file = "../MCMC_Runs/conjUVVsamp_minisimstudy_close_noSM_6_Dec_2023.rds")
# for(rep in 1:nreps){
# print(table(sim_results_UVV_no_sm[[rep]]$k))
# p1 = make_k_traceplot(k = sim_results_UVV_no_sm[[rep]]$k,
# group_assign = sim_results_UVV_no_sm[[rep]]$group_assign)
# print(p1)
# }with split-merge step
# sim_results_UVV_with_sm = parallel::mclapply(X = 1:10, mc.cores = use_cores,
# FUN = function(x){MVN_CRP_sampler_UVV(
# S = n_iter, seed = seeds[x], y = yreps[[x]]$y,
# alpha = 1, r = 10, g = 1, h = 50, nu = 2,
# fix_r = FALSE,
# mu0 = matrix(round((colMeans(matrix(unlist(yreps[[x]]$y), ncol = 2))),0), ncol = 1),
# lambda0 = diag(x = 15, nrow = 2),
# k_init = 1, diag_weights = FALSE,
# verbose = FALSE, split_merge = TRUE, sm_iter = 10
# )}
# )
#
# saveRDS(object = sim_results_UVV_with_sm, file = "../MCMC_Runs/conjUVVsamp_minisimstudy_close_withSM_6_Dec_2023.rds")This took too long…over 1.5 hours. Had to kill job for time. Need to reexamine and determine how to handle in future.